Scalar Field Dark Matter: head-on interaction between two structures 



o 
o 

(N 
> 

o 
m 

(N 
> 

oo 

O 

O ■ 

^ : 

Oh; 

6 ■ 
a: 



13 



Argelia Bernal 1 and F. Siddhartha Guzman 2 

1 Departamento de Fisica, Centra de Investigation y de Estudios Avanzados del IPN, AP 14-740,07000 Mexico D.F., Mexico. 
2 Instituto de Fisica y Matemdticas, Universidad Michoacana de San Nicolas de Hidalgo. Edificio C-3, 
Cd. Universitaria, C. P. 58040 Morelia, Michoacdn, Mexico. 
(Dated: February 5, 2008) 

In this manuscript we track the evolution of a system consisting of two self-gravitating virialized 
objects made of a scalar field in the newtonian limit. The Schrodinger-Poisson system contains 
a potential with self-interaction of the Gross-Pitaevskii type for Bose Condensates. Our results 
indicate that solitonic behavior is allowed in the scalar field dark matter model when the total 
energy of the system is positive, that is, the two blobs pass through each other as should happen 
for solitons; on the other hand, there is a true collision of the two blobs when the total energy is 
negative. 
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I. INTRODUCTION 

The most widely studied dark matter hypothesis con- 
sists in assuming that it is made of point-like cold parti- 
cles that are responsible for the formation of structure in 
the universe; among the most studied candidates nowa- 
days are the supersymmetric particles that would behave 
as a cold fluid made of particles. However, two problems 
associated to the point-like nature of dark matter are 
that the resulting gravitational collapse shows a central 
density profile that is not flat and on the other hand it 
predicts a non-observed amount of small structures. An 
alternative to ameliorate these two problems consists in 
assuming that the dark matter is made of an ultralight 
spinless particle, the so called Scalar Field Dark Matter 
Model (SFDM). In the cosmological frame, the analysis 
of such hypothesis indicated that the power of structures 
could be controlled through a parameter in the model, 
that is, the mass of the scalar field representing the spin- 
less particle [l|, H, H, 0] . Once the mass to of the boson 
is fixed the power spectrum suffers a cutt-off according 
to the mass of the smallest structure desired. An in- 
teresting assumption in such analysis is that the scalar 
field potential was a cosh-like potential, that behaved as 
an exponential at early times and as a free field case 
(quadratic potential) at late times, whose behavior was 
that of the usual cold dark matter model. Moreover, it 
was found that the SFDM enjoys the same advantages 
at cosmic scale as the standard lambda cold dark matter 
model. 

Because the SFDM requires the existence of a funda- 
mental scalar field for its reliability, it is natural to con- 
sider that this scenario fits very well within unification 
theory scenarios and braneworld models @. This by it- 
self is a good enough reason to consider the SFDM as an 
alternative powerful model, because it contains intrinsi- 
cally the spinless boson as dark matter particle. However, 
once at cosmic scales the model matches with observa- 
tions, it is necessary to study the predictions of the model 
at structure scales. In this sense there have been several 
results indicating that the model is good also at galactic 



scales and here we briefly summarize such results. 

At the early stages of the galactic dark matter model, 
fully general relativistic stationary solutions were pro- 
posed to explain the phenomena like the flatness of ro- 
tation curves, assuming the scalar field was real [f| 0], 
non-topological scalar field dark halos [H complex scalar 
fieldsQ and global monopoles [l(|. Later on, the as- 
sumptions relaxed to the newtonian limit of such general 
relativistic models. By the time, quintessential dark mat- 
ter halos [ill, El, El an d the fluid dark matter made of 
scalar fields w as p roposed as an alternative galactic dark 
matter modelpjj and the collapse of fuzzy dark matter 
made of a scalar field was analyzed in one dimension [TEl]. 
On the other hand, the assumption of time indepen- 
dence was also relaxed and scalar field dark matter halos 
were proposed to be gigantic Oscillatons, that is, time 
dependent fully relativistic scalar field solutions to the 
Einstein-Klein-Gordon system of equations 

Currently what appears to be the interesting case is 
that of the time dependent newtonian limit of the model, 
that is, the Schrodinger-Poisson (SP) system of equations 
would describe the model at local scales. In this direc- 
tion relevant results have been found, for instance: it was 
shown that when the evolution of a structure of galactic 
mass is followed after the turnaround point, it quickly 
virializes and tends to a stationary equilibrium solution 
of the SP system of equations, whereas one of the size 
of a supercluster would still be relaxing at the present 
time[l8|]; the condition is that the mass of the boson 
(to ~ 10 _23 eV) is the one that better cuts-off the power 
spectrum at galactic scales as shown in [l|, [H, HI- Thus, 
at the moment the pieces of the model seem to match 
both, at cosmic and at local scales. In fact, recently 
in [19( it was shown that the scalar field gravitational 
collapse tolerates the introduction of a self-interaction 
term in the potential, which makes the model to seem 
quite like a self-gravitating Bose-Condensate. In [2(1 we 
showed that spherically symmetric equilibrium solutions 
of the SP system are stable against non-spherical per- 
turbations, and moreover, such configurations played the 
role of late-time attractors for initially quite general ax- 
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isymmetric initial density profiles. 

What we present here is a step forward in the direc- 
tion of studying the evolution of scalar field structures. 
We perform numerical studies of scalar binary configura- 
tions, as a first step towards the making of a numerical 
code with no symmetries and for N-scalar objects. These 
studies would tell us about possible restrictions on self- 
interaction terms for the scalar field, and the way single 
configurations interact with each other. We restrict our- 
selves to the case of head-on interaction, which can be 
handled with a 2D code with axi-symmetry. We choose 
to write down the SP in cylindrical coordinates: 
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where tp — tp(x,z,t) and U = U(x,z,t) are the wave 
function and the gravitational potential respectively; x, z 
are the radial and axial cylindrical coordinates respec- 
tively. The third order term in Eq. (JTJ is related to a 
self-interacting term, in which A corresponds to the s- 
wave scattering length in the Gross-Pitaevskii approxi- 
mation for Bose Condensates [2l[. This term was shown 
to play the role of determining the compactness of the 
structure [ijj]. Equations |T][2]) use the units and scal- 
ing h = c = 1 with x — > mx, z — > mz, t — > mt and 
the wave function %p — > V AnGtp, where m is the mass 
of the ultralight boson. As a consequence of this change 
of units is that the mass of a system will be in units of 
[M] = Mpjm as found also for fully relativistic boson 
star solutions [n| [HJ ; this implies that the value m sets 
the physical lenght and time scale of the configuraions 
evolved. This mass, together with the scaling relations 
of the Schrodinger-Poisson system [T3, [2(| HH are the ba- 
sics for transforming back to physical units the system of 
interest (see an exmple below). 

The paper is organized as follows. In the next section 
we briefly describe the code used. In section Unl we con- 
struct initial data containing two spherically symmetric 
equilibrium configurations along the z axis. In section 
IIVI we show the results for the head-on interaction of the 
two structures. Finally in section [V] we draw some con- 
clusions. 



II. NUMERICAL METHODS 

The evolution. The most common numerical technique 
for time-integrating Eq. (TIJ is implicit with alternating 
direction splitting of the evolution operator [24j, [25| . The 
reason for this is that the evolution operator is unitary. 
Nevertheless, we used such method in [23], where no 
need for splitting the operator on the right hand side of 
Eq. (TTJ) was needed; in [l^, [2(| it was shown that explicit 
methods preserve also the number of particles and no 



significant difference in the results is found after using 
one method or the other. For the present case an explicit 
approximation of the full implicit method (in practice, a 
modified iterative Crank-Nicholson method 17]), with 
second order finite differencing to calculate the spatial 
derivatives is used. The reason to avoid using the 
implicit method is the difficulty to reduce the evolution 
operation to a tridiagonal system of equations when 
considering a non-zero A, which makes the Schrodinger 
equation a non-linear one, a situation not descussed in 

Poisson equation. Equation © is an elliptic equa- 
tion for U which we solve using the 2D five-point sten- 
cil for the derivatives and a successive over-relaxation 
(SOR) iterative algorithm with optimal acceleration pa- 
rameter (see e.g. [26| for details about SOR). In order to 
impose boundary conditions we made sure the bound- 
aries were far enough for the mass M = J \ip\ 2 d 3 x to be 
the same along the three faces of the domain and used 
the monopolar term of the gravitational field; that is, 
we used the value U = —M/r along the boundaries with 
r = \/ x 2 + z 2 . At the axis we demanded the gravitational 
potential to be symmetric with respect to the axis. 

We use a sponge in the outermost region of the do- 
main. The sponge is a concept used with success in the 
past when dealing with the Schrodinger equation (for de- 
tailed analyses see [H, HtJ)- This technique consists in 
adding up to the potential in the Schrodinger equation 
an imaginary part. The result is that in the region where 
this takes place there is a sink of particles, and therefore 
the density of probability approaching this region will be 
damped out, with which we get the effects of a physically 
open boundary. 

Basic testbeds of this code evolving single equilibrium 
configurations can be found in [20], where the results 
are also compared with previous studies with spherical 
symmetry and linear perturbation theory. 



III. INITIAL DATA 

Details about the construction of initial data for spheri- 
cally symmetric equilibrium configurations can be found 
in [19|, [2(| [23|]. Here we briefly mention the procedure 
used in our binary case. What we do is to superpose two 
spherically symmetric ground state equilibrium configu- 
rations upon the same 2D axially symmetric grid, whose 
construction is described as follows. In spherical symme- 
try equations (|1I2|) read 



id t ip 
d 2 (rU) 



1 

~27 

ripip 



d 2 r {riP) + UiP + K\iP\ 2 iP (3) 
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where r — v ' x 1 + z 2 . If a time dependence of the type 
if) — (pe lujt , regularity at the origin, and an isolation con- 
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dition <fi(r — > oo) = are assumed, the system becomes 
an eigenvalue problem for <f> with eigenvalue u>: 



d 2 r (r4>) 
d 2 (rU) 
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(5) 
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In order to solve these equations we discretize them and 
use a shooting method that bisects the value of u> so that 
the boundary conditions hold with certain desired accu- 
racy. The solutions constructed in this way are called 
equilibrium configurations. For each value of A it is pos- 
sible to construct the whole branch of equilibrium con- 
figurations as shown in [l9| . However there is an extra 
ingredient in these solutions: the number of nodes of the 
wave function. When the wave function is nodeless we 
say it belongs to the ground state; when the wave func- 
tion has a given number of nodes we can construct also 
the branches of solutions that -by analogy with the parti- 
cle in a box- are called excited spherical states. However, 
as shown in [19), 23] such excited states are unstable and 
decay into ground states; in fact ground states are stable 
and late-time attractors for quite arbitrary initial wave 
function profiles [H, H3| ■ Because the excited states are 
unstable and decay in a short time scale we collide in the 
present task only ground state configurations. 

Once we account with these ground state data: i) we 
interpolated the wave function of such configuration cen- 
tered at (0, Zo), ii) we place another of these equilibrium 
configuration at the point (0, — zq), iii) we choose zq so 
that the two configurations are far enough one from the 
other (see below) and iv) resolved Poisson equation @. 
Then we have initial data for two ground state equilib- 
rium configurations in our axially symmetric domain. 

Summarizing, we choose to solve the initial value prob- 
lem in spherical coordinates to make sure that we start 
the evolution with very accurate values, and we evolve in 
a 2D grid using cylindrical coordinates because we found 
them necessary and practical for the binary case we are 
interested in. 

Special warning is needed in Eq. {1} , because it is non- 
linear for the A ^ case, which indicates that the super- 
position of two wave functions is not allowed. Assume 
if>i and ip2 represent the solutions of the initial spheri- 
cally symmetric configurations that are to be superposed 
onto the 2D grid; the density of probability in ^ for the 
total wave function ip — ipi + -02 is \ipi + ip2\ 2 and un- 
less these states are orthogonal one cannot consider the 
naive superposition above is allowed. Thus we choose 
the distance between the configurations such that the in- 
terference (given by the scalar product of the two wave 
functions) < ipx,ip2 > is of the order of the precision of 
our calculations, say, in our case, the precision of the in- 
terpolation of the data into the 2D grid. Thus we can 
think of the system as one made of two adequately su- 
perposed equilibrium configurations we want to collide. 
An example of the interference term is shown in Fig. [TJ 

The superposition of configurations by itself would say 
little about whether or not two configurations collide. 



We add an extra ingredient to the system, that is, an 
initial head-on momentum to the initial scalar field balls. 
We simply generate different initial kinematical states by 
assigning new values to the wave functions of the equilib- 
rium configurations: ipi — > -0 1 e 4PiZ and tp2 — * 4>2e~ lPzZ - 
The resulting physical situation involves a considerable 
change in the value of the expectation value of kinetic 
energy in the system. 




FIG. 1: (Left) An example of two ground state configurations 
superposed along the z— axis with two different separations. 
(Right) The interference < ipi,ijj2 > is shown for two cases. 
Different separations are used and the two configurations have 
a central field value ip(0) = 1.0. The initial head-on momen- 
tum is p z — 3.0. 



Fortunately, in the present Newtonian low energy 
regime it is possible to estimate expectation values of 
physical quantities, a property difficult to pose in the 
fully general relativistic case. For instance, we account 
with the observables that allow one to monitor the evo- 
lution of the physical situation. Because we deal with 
a quantum mechanical system, we simply estimate the 
expectation values of the following interesting operators: 



K 
W 
I 



i)*W 2 il)d 6 x 



| / iP*Uif;d 3 x 



(7) 
(8) 
(9) 



which are the expectation values of the kinetic, gravita- 
tional and self- interaction energies. These quantities are 
quite important at determining the state of the system at 
any time during the evolution of the system. That is, the 
value of the total energy E = K+W+I indicates whether 
we account with a bounded system or not, and the very 
important virial theorem relation 2K + W + 31 = [28| , 
which is nearly satisfied when the system gets virialized 
and relaxed through whatever channels available, for in- 
stance, the emission of scalar field bursts, the so called 
gravitational cooling. 
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IV. HEAD-ON COLLISIONS AND SOLITONIC 
BEHAVIOR 



A. Equal mass case 



The first scenario one might think of is the collision of 
two equal mass ground state configurations. In Fig. [2] 
we show snapshots of the density profile along the z-axis 
for an initial configuration with p z = 3.0 and zq — 15 in 
the free field case (A = 0). What is found is that the 
two blobs move toward each other and eventually they 
lie upon each other, an interference pattern gets formed 
and after a while one blob moves toward the left and the 
other one toward the right. The first interpretation is 
that the initial data behave like solitons. Unfortunately 
we cannot be confident about the solitonic behavior be- 
cause the shape of the blobs gets deformed after the "col- 
lision" and an increase of the amplitude and shrink in the 
width arc manifest. After the distributions approach the 
boundary (located at z = ±30) the density of probabil- 
ity is absorbed by the sponge and its integral M drops 
to zero. At this point we are unable to track the evo- 
lution further in time and we ignore whether the blobs 
might return and collide again and repeat such process as 
many times as desired until there is energy released (e.g. 
through the emission of scalar field) and the encounters 
get damped allowing eventually a true collision. In Fig. 
[3] we show a zoom of the interference pattern at the time 
when superposition of the configurations around t ~ 5 
occurs. 

Of course, not all the initial configurations constructed 
present this behavior and we have found that a crite- 
rion to decide whether this behavior is allowed or not 
is the value of the total energy E = K + W + I. In 
Figs. 0] and [5] we show the total energy of different types 
of initial configurations. In Fig. 0] we present different 
situations for the free field case A = and two partic- 
ular cases: p z — 1.0 and p z = 0; in the first case the 
solitonic behavior is achieved and the total energy is al- 
ways positive and approaches zero because the density 
of probability has left the numerical domain; in the sec- 
ond case the total energy is always negative and at the 
end of the day what is found is that there is a single 
blob in the middle, indicating that the system is oscil- 
lating around a bounded object. We show snapshots of 
this behavior later when dealing with the more interest- 
ing unequal mass case. About the other cases in this 
plot p z = 0.75,0.71,0.7 we cannot decide whether they 
show solitonic behavior or not in the time we used to run 
our simulations and we can only observe that the density 
profiles are severely distorted by the collision. 

In Fig. [S]we show the same criterion for configurations 
with the self-interaction term (A = 0.2). Again, when the 
total energy is extremely positive or extremely negative 
we are able to decide whether the configurations collide 
or they trespass each other. 
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FIG. 2: Snapshots of the density of probability for interac- 
tion between two initial configurations with A = 0, p z — 3.0 
and zq = 15. The configuration shows solitonic behavior due 
to the fact that the total energy is always positive of the or- 
der of E ~ 32 until the blobs get absorbed by the sponge. 
The numerical domain used is x G [0,30], z £ [—30,30] with 
resolution Aa; = Az — 0.125. 



B. Unequal mass case 

The unequal mass case helps at deciding whether the 
configurations in the above examples truly trespass each 
other or bounce. In fact up to now it is not possible 
to say anything about this because of a few reasons: i) 
the expectation absolute value of the linear momentum 
along the head-on direction is equal for both half planes 
z > and z < 0, ii) the mass in each half plane is also 
the same in both semiplanes, iii) the expectation value 
of the linear momentum in the head-on direction is zero 
all the time. 

In the unequal mass case we have the advantage of 
being able to distinguish the half planes masses and linear 
momentum. In Fig. Owe show snapshots of the unequal 
mass case for A = 0.2 and initial parameters p z = 3.0 
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FIG. 3: The pattern of interference formed during the colli- 
sion of two configuration with A = 0, Zo = 15, p z =3.0 (the 
same case as in the previous figure). After this stage, the two 
blobs continue their way in the initial direction and behave 
as solitons (see Fig. [2j. They are not strict solitons because 
they slightly deform during the process. 

and Zq = 15 that show solitonic behavior with E > 
all the way. It can be seen clearly that the initial blobs 
are actually trespassing each other although they suffer 
a profile deformation until the wave function reaches the 
sponge region. 

In Fig. [7] we show the mass transfer from the z > to 
the z < half planes and viceversa; notice that the mass 
transfers from one side to the other in a very effective 
way. We also show the expectation value of the linear 
momentum along z in both half planes; what is found is 
that the momentum is also transfered from one side to 
the other. The perfect solitonic behavior would consist 
of having these two properties plus the unachieved one 
related to the preservation of the density profile. 

C. An example of collision 

As a final result we show what happens when an initial 
configuration presents a negative total energy. In Fig. [8] 
we show the density profile along z for a collision case 
corresponding to A = 0, zq — 10 and p z = 0, that is, 
only the pure gravitational force drives the dynamics of 
the binary. It can be observed that the blobs actually 
merge and remain sitting on a fixed point around the 
center of mass, the density tends to get stabilized, the 
virial relation starts oscillating aorund zero with smaller 
amplitude, the total energy starts stabilizing, so as the 
mass of the system. 

In order to ullistrate what our results mean in physical 
units we use the run in Fig. [5] to estimate the time-scale 
for the collision of binary equal mass head-on case. We 
start from the fact that the mass of ground state config- 
urations is M ~ 10 n M© and the mass of the boson is 
m = 10 _23 eV; the separation is 20 ~ 3.52kpc and the 
time the density peak is maximum is t co m S i on ~ 52.5 ~ 



FIG. 4: The total energy E = K + W + I for different ini- 
tial parameters, all of them with A = and zo = 10. The 
configurations with p z = 2.0, 1.0 start evolving with high val- 
ues of the total energy, and show solitonic behavior; the total 
energy tends to zero once the density of probability (the part- 
ciles) gets out of the numerical domain, which happens by the 
time indicated with the star for such run. The configuration 
with p z = remains with clearly negative total energy, so 
that the system is bound and the system collides, in the sense 
that there is no solitonic behavior and instead the two blobs 
get glued and remain like that. We are unable to conclude 
anything about the borderline cases. The stars at the end 
of borderline cases indicate the point at which we stoped the 
runs. The reason is that the lenghtscale of such cases is pretty 
much that of our numerical domain, and one expects the trun- 
ing points of the blobs to be at a distance of the order of the 
domain size. By the time indicated with a star, there is burst 
of particles, related more to the fact that the blobs are re- 
turning bach to the domain than to a burst of particles due 
to the relaxiation of a single blob or to the fact that the blobs 
are leaving completely the domain. 
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FIG. 5: The total energy E = K + W + I for different initial 
parameters, all of them with A = 0.2 and zo = 15. Once 
again, the configuration with a clearly positive total energy 
(p z = 1.0 and the unequal mass case with p z = 3.0) shows 
solitonic behavior and those with clearly negative total energy 
(p z < 0.5) show a merger. Once againg, we stopped the run of 
the borderline case with p z = 0.755, which physical properties 
change at about t ~ 50, and the numerical domain used does 
not suffice to determine its fate. 
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FIG. 6: Snapshots of the density of probability for the in- 
teraction between two unequal mass initial ground state con- 
figurations with A = 0.2, p z = 3.0 and zq = 15 (the first 
case in the previous figure). The configuration shows soli- 
tonic behavior due to the fact that the total energy is always 
positive of the order of E ~ 35 until the blobs get absorbed 
by the sponge. The total energy remains positive all the time, 
until the blobs are absorbed by the sponge. The numerical 
domain used is a; £ [0, 30] and z £ [—30, 30] with resolution 
Ax = Az = 0.125. 



8.3 x I0 6 yr; the maximum relative speed before the col- 
lision is v ~ 830fcm/s. 



V. CONCLUSIONS 

We have presented numerical solutions of the 
Schrodinger-Poisson system of equations which includes 
the non-linear term related to the self-interaction in the 
mean field Gross-Pitaevskii equation for Bose Conden- 
sates. In such case, the potential well is given by self- 
gravity of the density of probability of the system. The 
particular case we have studied corresponds to the inter- 
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FIG. 7: Top: we show the transfer of particles from the z < 
half plane to the z > one and viceversa. The mass in each 
semiplane is shown. Bottom: we show the transfer of mo- 
mentum between the two half planes, actually the expectation 
value calculated in each semiplane. The configuration evolved 
consists of two superposed ground state configurations with 
A = 0.2, and respective central field values ip(0, 0) = 1,0.7 
and masses M = 2.437, 1.9383. The initial linear momentum 
is p z — 3.0. The density of probability reaches the edges of the 
numerical domain and vanishes. The total energy is always 
positive. 



action between two ground state configurations (spheri- 
cal both of them) [13 ]. 

We found that the initial blobs show solitonic behavior 
of the initial configuration, but that also the two con- 
figurations may collide. The system ends up colliding 
whenever the total energy of the system E < and the 
solitonic behavior appears when E > 0. Unfortunately 
we can show this only for clearly non-zero values of the 
energy in each case and we cannot conclude anything 
about the borderline case, that is, when E ~ 0, because 
our simulations are unable to resolve the system for the 
quite long time needed and the spatial domain used. 

Within the scalar field dark matter paradigm, the two 
initial blobs would represent two virialized structures 
made of dark matter. What we have shown is that not all 
couples of configurations are allowed to have a collision, 
and that the total energy would indicate whether or not 
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a collision can occur. Our calculations also involve the 
presence of a self-interaction term in the scalar field, and 
therefore are within the Gross-Pitaevskii frame of Bosc 
Condensates, which this time are gravitating. 
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